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Abstract 

We study the nucleon matrix elements of the quark scalar-density operator using 
maximally twisted mass fermions with dynamical light {u,d), strange and charm 
degrees of freedom. We demonstrate that in this setup the nucleon matrix elements 
of the light and strange quark densities can be obtained with good statistical ac- 
curacy, while for the charm quark counterpart only a bound can be provided. The 
present calculation which is performed at only one value of the lattice spacing and 
pion mass serves as a benchmark for a future more systematic computation of the 
scalar quark content of the nucleon. 
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1 Introduction 



A number of experiments (see for instance refs. [IHS]) have been designed to 
investigate the nature of Dark Matter (DM) by detecting and/or measuring 
the recoil energy of nuclei hit by a hypothetical DM particle. One very pop- 
ular example of such DM candidates is a weakly interacting massive particle 
(WIMP) like the ones that are predicted in a large class of models [6]. However, 
definite evidences for a direct detection of WIMPs have not been observed up 
till now. Nonetheless, the various ongoing experiments provide rather severe 
constraints for the parameters of many DM models. 

A possible scenario for the detection of a WIMP type of DM particles relies 
on the idea that the WIMP - due to its assumed large mass - produces a 
Higgs boson which in turn couples to the various quark flavour scalar density 
operators taken between nucleon states, as depicted in Fig. [TJ In fact, at zero 
momentum transfer, the cross section for spin independent (SI) elastic WIMP- 
nucleon (x^) scattering reads [7] 



C^SI.xN 



T.G,fH)fTf with f^^ = ^{N\qjqj\N) . (1) 
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The functions Gg^ depend on several parameters of the particular model Cgi^^ 
is computed in, among which the WIMP mass, rriy.. A precise expression of 
the cross section in the constrained minimal supersymmetric Standard Model 
(CMSSM) can be inferred from ref. [7]. The expression above for ctsi^^n is in 
particular valid in the SU(2) isospin limit for u and d quarks, an approximation 
that is always understood in the following. What we want to compute is the 
magnitude of the dimensionless and renormalization group invariant (RGI) 
coupling J'tj:, which depends on the mass m^^. of the quark of flavour / and 
the nucleon mass, rriN. For short in the following we will refer to the matrix 
element {N\qfqf\N) as the "scalar /-quark content" of the nucleon. 

As one can see from Eq. ([T]), the cross section, ctsi^^^n depends quadratically on 
fTf , and it is thus very sensitive to the size of the scalar content contributions 
of different flavours. Already a 0(10%) variation of /t^ can lead to significant 
changes in 0"si,xN- It is therefore necessary to compute accurately and with 
controlled error the hadronic matrix elements {N\qfqf\N) . The implications 
of the hadronic uncertainty in comparing models with present data of WIMP 
direct detection has been for instance presented recently in [S]. 

One way to calculate for various flavours the scalar quark content of the nu- 
cleon is provided by chiral perturbation theory (xPT). The results of such a 
calculation are usually parametrized in terms of ct^at and ctq defined by the 
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Fig. 1. The Higgs-boson exchange contribution to the WIMP-Nucleon low energy 
scattering process. 

formulae 

a^N = mi{N\uu + dd\N) and ao = mi{N\uu + dd - 2ss\N) , (2) 

where rrii is the common mass of the u and d hght quarks. A low energy 
theorem relates aj^N to the pion nucleon scattering amplitude extrapolated to 
the Cheng-Dashen point [9]. The functional form of the extrapolation formula 
can be established using dispersion relations and in this way a-^N has been 
found to be ct^at = 79(7) MeV [TU]. Note that a recent result which uses Lorentz 
covariant baryon chiral perturbation theory gives a.^N = 59(7) MeV [TT]. The 
non-singlet matrix element ctq can be obtained by looking at the pattern of 
the SU{3) symmetry breaking visible in the spectrum of the baryon octet. A 
determination using xPT obtained in |T2] has given ctq = 36(7)MeV. 

A direct measure of the magnitude of the strange quark content of the nucleon 
relative to the light quark content is represented by the ratio 



The ?/Ar-parameter can be related to aT,N and ao by t/at = 1 — cto/ct^at and 
one finds, using the previously quoted numbers, i/n = 0.44(13) [13]. This 
result leads to a surprisingly large strange quark content of the nucleon, and 
consequently to a large strange quark contribution to the o"si^i^n cross section 
in Eq. ([T]). On the other hand, the quoted uncertainty on is rather large, 
i.e. of the order of 30%. As pointed out in ref. [2] and discussed above, it is 
quite important to provide a precise value for the strange quark content of the 
nucleon in order to be able to interpret the ongoing and planned experimental 
searches of WIMPs. 

Lattice QCD can provide a determination of these nucleon matrix elements 
from first principles. The difficulty involved in these computations has limited 



Vn = 



2{N\ss\N) 



(3) 



{N\uu + dd\N) ■ 
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for a long time the possibility to isolate a physical signal for these quantities. 
The aim of this paper is to show that it is indeed possible to accurately 
compute uiq and test whether the strange quark content of the nucleon is as 
large as indicated by xPT. To this end, we will use Wilson lattice QCD with 
maximally twisted mass fermions since this framework is particularly suited 
for such a calculation, as we will discuss below. 

On the lattice, two approaches have been followed for the calculation of the 
light and strange quark content of the nucleons. The first is based on the 
Feynman-Hellman theorem [T3] that relates the nucleon scalar matrix element 
to the dependence of the nucleon mass on the quark masses via the equations 



cr^^ = mi{N\uu + dd\N) = mi , (4) 

omi 

= ms{N\ss\N) = rus^^ , (5) 

where we denote by the strange quark mass and the derivative are intended 
to be evaluated at the quark mass values that correspond to the physical pion 
and kaon masses. This approach is often referred to as the "spectrum method" . 
The other approach, the so-called "direct method" , consists in evaluating di- 
rectly the matrix elements appearing in Eqs. (jl]) and ([5]). 

Both approaches are numerically very challenging. In the spectrum method, a 
number of simulations at different values of the strange quark mass are needed 
in order to evaluate the derivative. In the direct method, the computation of 
the disconnected diagrams is required, which is a highly demanding task since 
they often show quite a bad signal to noise ratio, thus requiring very high 
statistics. In addition, as pointed out in ref. [I6], when lattice discretizations 
that break chiral symmetry are used, a mixing between the bare light and 
strange scalar quark density matrix elements occurs under renormalization. 
This mixing is not present for chiral invariant (e.g. overlap) fermions |17j . 
It can also be avoided up to O(a^) effects when - like we do in this work - 
maximally twisted mass fermions are used in a mixed action setup as explained 
below. 



There exist already a number of lattice QCD computations of the strange 
quark content of the nucleon, see the works of refs. [T6l[T8l - [28] and the review 
paper [13]. Although these calculations indicate that the strange quark con- 
tent of the nucleon is smaller than suggested from xPT, the quoted results 
are affected by large statistical errors, making it difficult to reach definite 
conclusions. 

In this paper we present a method which allows to compute the strange quark 
content of the nucleon with small statistical errors. By means of a benchmark 
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calculation at only one value of the lattice spacing, volume and quark mass we 
demonstrate that it is indeed possible to achieve for the parameters /t^ and Un 
(see Eqs. ([T]) and (|3])) a signal to noise ratio significantly (more than 4 standard 
deviations) different from zero. It remains of course open the question of the 
size of systematic effects which we plan to address in the future by including 
in the computation more lattice spacings, volumes and quark masses. 

As we said, a main ingredient that enables us to obtain accurate values for the 
strange quark content of the nucleon is the use of Wilson lattice QCD with 
(maximally) twisted mass fermions for our simulations. Besides its property of 
automatic 0(a)-improvement, it offers the advantage that special techniques 
for computing disconnected diagrams (see below) can be employed, which 
significantly reduce the size of the numerical noise typically affecting the com- 
putation of such diagrams. In addition, with our choice of twisting for valence 
fermions the operator matrix elements relevant for the various flavour contents 
of the nucleon turn out to be all multiplicatively renormalizable. Although the 
same property is valid in chiral invariant (e.g. overlap) lattice formulations, 
twisted mass fermions are computationally much less demanding, enabling us 
to work with a large number of gauge configurations and at fairly big volumes 
and small lattice spacing. 

As a last point, we want to mention that for this computation we employ 
gauge configurations generated with Nf = 2 + 1 + 1 dynamical quarks, in which 
besides a mass degenerate u and d light doublet also a mass non- degenerate 
strange s and charm c pair is present in the sea. This will allow us for the first 
time to also study the charm quark content of the nucleon. 



2 Computational methods 

2.1 Lattice action 

The lattice action used in our simulations includes as dynamical degrees of 
freedom, besides the gluon field, a mass-degenerate light up and down quark 
doublet as well as a strange-charm quark pair, a situation which we refer to as 
the Nf = 2 + 1 + 1 setup. While in the pure gauge sector we use the Iwasaki 
action [29], for the fermion part twisted mass fermions are used. In particular, 
concerning sea quarks, we make use of the formulation of refs. [SUIEI] for the 
light mass degenerate u-d sector, while the action introduced in refs. [521153] is 
employed for the mass non-degenerate c-s sector. The quark mass parameters 
of the heavy flavour pair have been tuned so that in the unitary lattice setup 
the Kaon and D-meson masses, take (approximatively) their experimental 
values. More information about this scheme and further simulation details 
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can be found in ref. |34] . 



In order to fix the notation, we give here the exphcit form of the twisted mass 
action for a doublet of degenerate quarks : 



1 — 7 1 + 7 il 

-Xfix) E ^Uf,{x) — ^Xfix + a/i) + Ulix - afi) ^ ^ Xf{x - a/i)J \ . (6) 

Here /i/ > denotes the bare twisted mass. The hopping parameter is an 
alias for the bare standard mass mo/ = ((2k/)-i -4)/a. Eq.E] defines d[{^[U] 
the flavour degenerate Wilson twisted mass operator. When is tuned to its 
critical value, the lattice QCD formulation known as maximally twisted 
mass fermions is achieved, which guarantees 0(a) improvement of physical 
observables. The value of Kcr for all valence flavours is taken to be the same 
as in the sea quark sector 



For further needs we also introduce the operators -D/,± denoting the upper 
and lower flavour components of d[1^\U\, referred to as the Osterwalder-Seiler 
Dirac operator : 

DfAU]=t^ {i^D«[U]}, (7) 

where tr denotes the trace in flavour space. Df^±\U] then is the Dirac operator 
of an Osterwalder-Seiler lattice quark with mass ±/i/. 

When we discuss below the 2-point and 3-point correlation functions necessary 
for this work, we will use the so-called physical basis of quark fields denoted 
a.s ipf. The physical field basis is related to the twisted quark field basis, Xfi 
by the following field rotatioiJXl 



e 



ipf = ^^'■^ Xf and 'ipf = Xf 

where the twist angle uj = 7r/2 at maximal twist. In the following, ijjf with 
index / = /,s,c will denote quark field doublets of light (/), strange (s) or 
charm (c) quarks depending on the mass fif chosen in the valence sector. 
Since ipj will always refer to the physical basis we will denote with u and d 
the two components of ipi. Staying close to the notation of Eq. ([7]) we will 
denote with s± (resp. c±) the two components oi ips (resp. ipc)- 



^ At maximal twist the quark fields 'il)f,il!f are said to be in the "physical" quark ba- 
sis if the quark mass term in the Lagrangian appears in the canonical form ip ^fijipj. 
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In order to have a consistent mixed action setup, the values of the bare quark 
mass parameters Hs and fic in Eq- © have been tuned such that the Kaon 
and D-meson masses of the unitary setup [M] are matched. 

In particular the matrix element entering (Tttat will be calculated in such a way 
that (Euclidean) unitarity is preserved at finite lattice spacing, while those 
of interest for the strange and charm content of the nucleon are evaluated in 
a mixed action setup where unitarity violations represent mere O(a^) arte- 
facts [33]. Such cutoff effects are expected to be numerically small in line 
with the findings from previous mixed action studies carried out on ETMC 
Nf = 2 [32] and Nf = 2 + 1 + 1 [SB] gauge ensembles. 

2.2 Nucleon scalar matrix elements 

The nucleon two-point function is defined in the physical quark basis by 

= {v^{Jn{x)J^{x,,,))] , (9) 

X 

where (...) denotes field average, F^*" = (1 ± 7o)/2 are the parity projectors, 
Xsrc = (4rc5 s^src) is the spacc-time location of the source and r = t — tgrc stands 
for the source-sink separation. The subscript refers to the proton or to the 
neutron state for which the interpolating fields are given by the formulae 

Jp = e'''''(u''-^C^^(i')u' and Jn = e"^' (d"^^ C^i^u^) (f , (10) 

where C is the charge conjugation matrix. Note that due to translational invari- 
ance C^_2pts(''") does not depend on the spatial source location, Xsrc, which we 
can thus choose freely. Let us also recall that, since here we work in the SU (2) 
isospin limit approximation, an exact symmetry of the action (i.e. Vx{u d) 
where V is parity) leads to the relation C^2pt(''") = C*^2pt(''") [SZ]- 

The three-point functions of interest in this paper are defined by 

CUv^{t,^o^) = E tr {T^{JNix)Of{x,,)Mx,,,))} , (11) 

where Of for f = I, s,c denotes an appropriate (see below Eq. (|T3ll ) lattice 
regularization of the light, strange or charm quark scalar density and Top = 
^op — ^src is the operator-to-source time separation. Since we are considering 
an operator with a non vanishing vacuum expectation value, we also need to 
introduce the corresponding vacuum subtracted correlator 

Cpi7\T, Top) = Cpi,{T, Top) - C^,2pt(r) J:(0.M) ■ (12) 
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To be specific, the operator Of will be given for our case by 

1 _ _ 1 _ _ 

Oi = uu + dd, = 2 (■^+■5+ + and = 2 + c_c_) , (13) 

depending on the quark flavour of interest. Using these operators, we shall 
obtain the multiplicatively renormalizable, 0(a) improved matrix elements 
relevant for this paper. 

For each flavour /, the bare scalar matrix element at zero momentum trans- 
fer {N{p)\Of{0)\N{p)) can be written in terms of an effective coupling bare 
constant gsj and the nucleon spinor un, in the form 

(iV(p)|O/(0)|iV(p)) = gsjUN{p)uN{p) . (14) 

Using the two- and three-point correlators of Eqs. and (1 11 1) , we build for 
f = l,s,c the ratio 

Rf{r, Top) ^ "^^^g =gs,f + O(e-^^^l-^l) + O(e-^^^l--^l) , (15) 

where AM is the mass gap between the lowest nucleon state and the first 
excited state with the same quantum numbers. One can thus extract from the 
asymptotic time behaviour of the various -Rj(r, Top) the bare effective scalar 
couplings Qsj, which are in turn simply related to the nucleon sigma terms 
of interest. For instance, at maximal twist the lattice regulated versions of 
Eqs. (jl]) and ([5]) will read 

(T^N = ms, , cr^^t = ^,gs^ . (16) 

The systematic errors 0(e~^^''^''^~'^°p') and 0(e~^*^''^°p') originating from the 
finiteness of the time separations r — Top and Top will be neglected in this work. 
However they can have a non- negligible impact on the evaluation of nucleon 
matrix elements, as shown for instance in [38]. We are therefore planning to 
address this problem in a forthcoming publication. 

2. 3 Lattice discretization and evaluation of correlators 

The main aim of the paper is to study whether the improved methods to com- 
pute disconnected diagrams as applicable for twisted mass fermions will indeed 
lead to a calculation of the quark contents of the nucleon with significantly re- 
duced errors compared to earlier works. The analysis performed in this work 
concentrates therefore on one ensemble of a 32^ x 64 lattice volume with a 
lattice spacing of a = 0.0779(4) fm (/3 = 1.95) where the error quoted is only 
statistical [55], and a pion mass of approximately 390 MeV (a/x/ = 0.0055). 
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In order to improve the overlap between the ground state and the interpo- 
lating fields we use Gaussian smearing of the quark fields appearing in the 
interpolating fields. We also use APE smearing of the gauge links involved in 
the Gaussian smearing, following the same strategy as in 



In the twisted basis the scalar operators read 

Of = iXflb-r^Xf, where f = l,s,c , (17) 

and are hence given by the pseudo scalar density. There is a one-to-one cor- 
respondence between the bare operators Of and the bare operators Of intro- 
duced in Eq. f[T^ which is given at maximal twist by: 

r if / = / 

Of = zXfl5r'xf = i^fi^f=\20s iif = s. (18) 

[20, if/ = c 

While the two-point nucleon correlators of Eq. (Q give only rise to quark- 
connected Wick contractions, in general the three-point functions of Eq. f ill I) 
yield both quark-connected (illustrated in Fig. [2^) and quark-disconnected 
(illustrated in Fig. [2]d) contributions. In the following we will refer to them 
simply as to connected and disconnected fermionic Wick contractions (or di- 
agrams) and shall write 

^Ar,'3pt(''"' ''op) = ^Af'sptl''"' ''op) + ^Ar'sptl''"' ''op) (19) 

with C^'gpt (resp. ^^^^spt) corresponding to the connected (resp. disconnected) 
quark diagrams, defined as 



<^^it(r,rop) = E {^^{[jN{x)Of{x,^)Mx,,,)])} , (20) 
^^^'ipt(^,^op) = E {t^{[Jn{x)Mx,,,)] [Oy(xop)])} , (21) 



where the symbol [...] is a shorthand for all the connected fermionic Wick 
contractions. In particular, the contribution of the disconnected fermion loop 
to ^^^'3pt on a given gauge configuration U in our setup reads 



[OfM] = -iWf E tr I 75 ( -(j;^ 



(22) 

where, in view of Eq. (fT8|) . wi = 1, Ws = Wc = 1/2 and d[1^^\U] are the 
Osterwalder-Seiler Dirac operators defined in Eq. ([7]). 
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For the strange and charm content of the nucleon only disconnected dia- 
grams contribute to the three-point correlator, while for the light quark con- 
tent both kinds of fermionic diagrams matter. The connected contributions to 
{N{p)\Oi\N{p)) have been evaluated using standard techniques for three-point 
functions ("sequential inversions through the sink"). In this method one needs 
to fix the sink-to-source separation r = t — tsrc and we choose, as in ref . |10] , 
r = 12a corresponding in physical units to a separation of r ^ 1 fm. 

Since, using discrete symmetries and anti-periodic boundary conditions in the 
time direction for the quark fields, one finds 



where T denotes the lattice time extent, in order to increase the signal over 
noise ratio we have averaged contributions related by the symmetry rela- 
tions (l23l) and (l24l) . In addition, we have carried out Dirac matrix inversions 
at a number (denoted by A'^^rc in the following) of randomly chosen source 
points per gauge configuration, with the goal of better exploiting the gauge 
field information contained in each configuration. 



Fig. 2. Connected (left) and the disconnected (right) graphs arising from the Wick 
contractions of the 3-point function. 

An important issue to be discussed is renormalization of the correlators in- 
troduced in sect. 12.31 The technical arguments are given in Appendix A. We 
summarize here the conclusions. After having subtracted the mixing with the 
identity in the correlation function (see Eq. f lT^ ). the operator Oi, Og and Oc 
do not mix among each other. Since also the bare quark mass renormalizes 
multiplicatively with a renormalization constant that is precisely the inverse of 
the one occurring in the renormalization of Oj, the lattice quantities fi/Qsj , for 
/ = l,s,c, yield 0(a) improved renormalization group invariant (RGI) sigma 




(23) 
(24) 




(a) 



(b) 
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terms (see Eq. (fT6i) ). The renormalization pattern is thus as straightforward 
as for chirally invariant overlap fermions. 



2.J^ Numerical estimate of disconnected loops 



Let us shortly sketch the variance reduction method for the evaluation of 



tr < r — rr—, 7-r-\ / , T = some Dirac matrix , 

with twisted mass fermions introduced in . The method has already been 

applied to study the rj' meson in |13]. It relies on the fact that the difference 
between the twisted basis Dirac matrices d[{^'^^ and d[1^~^ is proportional to 
the identity , 

A^^^^ - Dit^ = 2z/in5 (25) 

implying that 



PfU^^) r,U^+) r,U,+) V^t"! ) p,{f,-) ^"'"^ p,{f,+) T^if,-) ■ 

J-^tm J^tm ^tm J-^tm J-^tm J-^tm 

(26) 

For the practical calculation, we introduce a set S of A'^ independent random 
volume sources, {^[i], . . . , ^[^j, . . . , ^[at^]}, satisfying 



lim 

Nc^oo 



eM(a:)*eWy)]^ = M^^' (27) 



where i = 1, 12 refers to the spin and color indices of the source and [■ ■ ■]e 
denotes the average over the A'^ noise sources in S . 

Multiplying Eq. (12^ by a F matrix and the noisy sources ^*r]{y) and ^ir]{x) 
and taking the trace over spin and color indices we get 



tr <;f(^-^] \+0(Nr'/' 



'tm -'-^tm / {x,x) 

(28) 

where 

<P[r] = {l/D[t''%r] and =^t_j(l/A^^+))t = 4]75(VA^m~V5. (29) 

For the generation of the random sources we have used a Z2 noise taking all 
field components randomly from the set {1,-1}. We note that in the case of 
F = 75, the quantity in Eq. (!28|) . after summation over x = Xop, provides an 
unbiased estimator of the disconnected fermion loops of Eq. 
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2. 5 Performance of the method for disconnected loops 

As a first step, we performed a study to determine tlie optimal number, A^^, of 
stocliastic volume sources to be used for evaluating the disconnected diagrams 
of Eq. fl22|) . For a given number Ncoai of gauge field configurations increasing 
beyond some value will not improve the signal over noise ratio (SNR) since 
the noise induced by the fluctuations of the gauge fields will eventually become 
dominating. 

For the present test we have used 677 gauge field configurations and chose 
fixed time separations, r = 12a and Top = 6a. Fig. [3] shows the SNR as a func- 
tion of the number of stochastic source samples employed to evaluate the 
disconnected loops of Eq. (122]) . As indicated in the figure, for the computation 
of the two point function we have used Agrc = 1, 2, 3, 4 randomly chosen point 
sources per gauge configuration. 

The figure demonstrates that the signal over noise ratio reaches a plateau for 
A^ > 7, meaning that for larger values of A^c the error is dominated by the 
gauge field fluctuations. The finding that for A^ > 7 a plateau of the SNR 
is reached holds true for all values of A"src we have used. In addition we have 
checked that the above conclusion remains essentially valid in the quark mass 
range we intend to explore. 

Fig. in] also demonstrates that the SNR increases when more source points 
per gauge field configuration are used. When we change the number of source 
points from 1 source per configuration to 4, we find a decrease of the error 
by a factor of approximately ~ 1.6. Although this does not correspond to the 
optimal factor of 2, using A'src- values moderately larger than 1 turns out to 
be a convenient and economic way to increase the signal over noise ratio. In 
the final analysis, we will always use A^ = 12 and A'src = 4. 

In Fig. |4]we compare the efficiency of the method discussed in sect. 12.4] which 
is based on the peculiar property, see Eq. (125]1 . of twisted mass fermions, with 
another noise reduction technique relying on the hopping parameter expansion 
of the Dirac operator. This latter technique is not restricted to twisted mass 
lattice QCD and has been introduced in P5]B6]. We refer the interested reader 
to the appendix B of [13] for an implementation in the case of twisted mass 
fermions. As can be seen from the figure the twisted mass specific variance 
reduction technique improves the signal over noise ratio by a factor ~ 3. 
Performing a simple extrapolation in the number of gauge field configurations 
we estimate that with the hopping parameter expansion technique (9(10000) 
configurations would be needed to reach a result 5cr away from zero, while 
only (9(1000) configurations are necessary to obtain the same accuracy with 
our twisted mass specific technique. 
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Fig. 3. Signal to noise ratio (SNR) of the quantity i?s(T, Top), see Eq. (fTSjl . for fixed 
values of r/o = 12 and Top /a = 6 as a function of and for different values of Ng^c- 
The bare strange valence quark mass is afis = 0.016. The number of configurations 
used is 677. 



3 Results 



3.1 The Pion-Nucleon a -term, a^^N 



We first concentrate on the determination of a-j^^ defined in Eq. ([2]). Since 
for this quantity only the up and down quarks come into play, we work in a 
fully unitary setup, where valence and sea quarks are regularized in the same 
way. In the following we will denote by -Rconn. (resp. -Rdisc.) the contribution 
of (5^;4t (resp. I^^'jpt), see Eqs. ([201 EI]) , to the ratio Ri defined in Eq. ffTSj) . 
In Fig. Owe show our results obtained for Ri, -Rconn. and -Rdisc. as functions 
of Top = top — tsrc for a fixed sink to source separation, t = t — tsrc = 12a. 
-Rdisc. has been computed using measurements over 842 configurations with 
= 12 randomly chosen volume sources. -Rconn. has been computed using 
510 configurations by employing the fixed sink method. 

The connected part, -Rconn., denoted by the black filled circles, shows a pro- 
nounced time dependence indicating the contribution of excited states. In this 
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Fig. 4. Signal to noise ratio (SNR) of the quantity Rs{l2a,6a), see Eq. (fT5|) . for 
= 12 and Nsrc = 4 as a function of the number of gauge field configurations 
-^conf 1 for the variance noise reduction technique used in this paper and the hopping 
parameter expansion technique. The dashed curves are drawn to guide the eye. The 
bare strange valence quark mass is the same as in Fig. [3l 

work, we do not attempt to quantify the size of this systematic effect since 
our goal here is more to investigate whether statistically significant values for 
the scalar quark contents of the nucleon can be obtained. The disconnected 
part, -Rdisc.5 denoted by blue triangles in Fig. [5], clearly corresponds to a small 
contribution compared to the connected part -Rconn., of the order of ~ 10% of 
the full ratio, Ri, represented by the red diamonds in the figure. 

In Fig. [6] (a zoom of Fig. [5]) we show only the disconnected contribution. Note 
the change in the scale on the vertical axis. It is encouraging that, by employ- 
ing the techniques described above, we can indeed obtain a non-zero signal 
at a ~ 4(7 level. In order to determine the "plateau" values of -Rconn. , -Rdisc. 
and Ri, we performed several fits to a constant through the data varying the 
fit interval. The results are summarized in Table [H We find that the discon- 
nected contribution is about ~ 8% of the connected one. Nevertheless since 
the error on the connected contributions is smaller than the value of -Rdisc, 
the disconnected contribution cannot be neglected when computing the ratio 
Ri. We finally remark that all the statistical errors in this work are computed 
using the bootstrap method [77] . 

An estimate of the systematic error on ct^tv can be given on the basis of the 
spread of the results one gets by varying the time interval [rop^ , Top^] over which 
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5 10 15 



Fig. 5. Plot of the contributions i?disc. (blue triangle), i?conn. (black circles) and of 
their sum, Ri (red diamonds) as function of Top at r = 12a for a^u; = 0.0055 and 
f3 = 1.95. 





-Rdisc. (842 meas.) 


-Rconn. 


(510 meas.) 


Ri 


b'op^ 1 Top2 ] 


fit 


X^/ndof 


CL 


fit 


X^/ndof 


CL 


fit 


X^/ndof 


CL 


[2, 10] 


0.66(16) 


4/8 


0.85 


9.6(5) 


12/8 


0.14 


10.3(5) 


15/8 


0.05 


[3,9] 


0.72(17) 


1/6 


0.98 


9.9(6) 


3.2/6 


0.78 


10.6(6) 


4/6 


0.65 


[4,8] 


0.76(18) 


0.26/4 


0.99 


10.2(6) 


0.6/4 


0.97 


10.9(6) 


0.75/4 


0.94 


[5,7] 


0.79(18) 


0.02/2 


0.99 


10.3(6) 


0.04/2 


0.98 


11.1(6) 


0.05/2 


0.98 



Table 1 



Plateau values for the ratio i?disc. , -Rconn. and Ri relevant for the extraction of a-^N 
for different time intervals [Top^jTop^]- We also include the by degrees of freedom 
(x^/ndof) and the confidence level (CL). 

the plateau is taken, as displayed in Tabled] For this purpose we construct the 
distribution of all fit results, weighted by their confidence level, and take the 
variance of this distribution as our estimate of the systematic errors. Using 
this procedure we find 

a^N{mps ^ 390 MeV) = 151(8)(4) MeV, (30) 

where the errors correspond to statistical and systematic uncertainties, re- 
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Fig. 6. Zoom of Fig. [5] only showing i?disc. versus Top at r = 12a. 

spectively. Note that the dominant contribution to the systematic error comes 
from the connected part of the ratio Ri. Also, the value obtained here cor- 
responds to only one pion mass of about 390 MeV and an extrapolation to 
the physical value of the pion mass will be finally needed. Notice that chiral 
perturbation theory predicts that ct^tv vanishes in the zero quark mass limit. 
While the calculation of ct^tv at several quark masses and lattice spacings is 
beyond the scope of this paper, we remark that simulations in this directions 
are under way. 



3.2 Strange content of the nudeon 



For each value of the valence quark mass, one can define the quantity: 

fTf = ^(mOflN) = f = l,s,c, (31) 

where /x/ is the bare quark mass and (A^|0/|A^) the bare matrix element 
corresponding to the quark flavour /. As argued in Appendix A, /^^ is a RGI 
quantity in our mixed action setup. For the ensemble used in this study, the 
nucleon mass am^v has been determined in [H] and is am^r = 0.510(7). In 
Fig. [7| we show -Raise, versus Top for a quark mass of a/^s = 0.016. We see that 



16 



the ratio is ~ 5a away from zero in the middle of the "plateau" . 
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Fig. 7. Plot of -Rdisc, i-e. the bare ratio of Eq. (fTSj) . versus Top at r = 12a in the 
strange quark mass regime {afig = 0.016). 

We have performed several fits varying the time interval to extract a "plateau" 
value. Results for /rs,o"o and i/n- are summarized in Table [21 While yiy does 
not depend strongly on the fit window and seems thus to be, within our 
accuracy, free of excited state contaminations, Jt^ and ctq are affected by an 
excited state contamination in a similar way as we observed for a.,rN- As we 
already remarked before, a quantitative evaluation of this systematic effect 
goes beyond the goal of this paper and will be addressed in the future. Our 
present results at a pion mass about 390 MeV for f^^, do and i/n are 



/t. =0.014(5)(1), (32) 
ao = 137(7) (4) MeV, (33) 
1/^ = 0.082(16)(2), (34) 

where the first number in parenthesis represents the statistical error and the 
second the systematic uncertainty. In order to estimate the magnitude of sys- 
tematic effects the same strategy as in the case of a-^N has been employed. The 
statistical error on /t^ also includes the error on the nucleon mass determina- 
tion. Note that the value of t/N quoted is obtained directly from the ratio of 
three point functions C'^,3p"^(r, Top) and Cp^''^p^^{T, Top) and agrees within error 
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with the value of un estimated from ao and cTttat {vn = 1 — o'q/o'-kN ~ 0.092). 
We stress again that our measurement of the strange content of the nucleon 
leads to a value of un which is different from zero at a Sex level. It is inter- 
esting to compare /^^ to the value, /t^, one gets in the light quark sector, for 
which we obtain a significantly larger value, namely fxi = 0.117(12)(3) (at 
mps ~ 390 MeV). 







(Jo (MeV) 


yAT-parameter 




fit 


X^/ndof 


CL 


fit 


X^/ndof 


CL 


fit 


X^/ndof 


CL 


[2,10] 


0.013(2) 
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11/8 


0.18 


0.08(2) 


0.7/8 


0.99 


[3,9] 


0.014(2) 


0.8/6 


0.99 


135(8) 


3/6 


0.79 


0.08(2) 


0.2/6 


0.99 


[4,8] 


0.014(3) 


0.2/4 


0.99 


138(8) 


0.5/4 


0.97 


0.08(2) 


0.03/4 


0.99 


[5,7] 


0.015(3) 


0.03/2 


0.99 


140(8) 


0.01/2 


0.98 


0.09(2) 


0.01/2 


0.99 



Table 2 



Results of the fits to /ts,o"o and at a/is = 0.016. The notations are the same as 
in Table [TJ 

3.3 Charm quark content 

Following the same strategy as in sect. 13. 2[ we have carried out the first study 
of the charm quark content of the nucleon. This is possible because we have at 
our disposal Nj = 2 + 1 + 1 simulations with a fully djTiamical charm quark 
degree of freedom. 

We show in Fig. |8]the dependence of Rc on Top (at r = 12a), using exactly 
the same statistics as in the light and in strange quark sectors. Unfortunately, 
for the charm quark content no hint of a plateau is visible. Signal and noise 
have equal order of magnitude and our results are compatible with zero. For 
comparison we also show the results for the strange quark content obtained 
in the previous section as a grey band. From our data we can only establish 
the inequality 

l(iV|a|iV)|<|(iV|Os|iV)|. (35) 



4 Conclusion 

In this work we have performed a benchmark calculation of the scalar quark 
contents of the nucleon by directly computing the matrix elements (A^|Oj|A^) 
for / = /, s, c. Extending the calculation to strange and the charm quark 
flavours became possible owing to Nf = 2 + 1 + 1 dynamical simulations 
recently carried out by the ETM Collaboration [53]. Our calculations were 
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Fig. 8. The ratio Rc of Eq. (llSp for fixed r = 12a versus Top in the charm quark 
mass regime (a//c = 0.2). For comparison, we display as a grey band the plateau 
value obtained in the strange quark case, see Fig. [71 

performed at one value of the lattice spacing (a ~ 0.078 fm) and for fixed 
values of the pion and Kaon mass (mpg ^ 390 MeV and rriK ~ 580 MeV, 
respectively) . 

In evaluating these nucleon matrix elements, maximally twisted mass fermions 
are very helpful in two respects. The first is that the twisted mass fermion 
regularization provides a framework where it is possible to efficiently evaluate 
quark-disconnected diagrams. The second is that a consistent lattice frame- 
work can be set up where the matrix elements of interest are multiplicatively 
renormalizable and at the same time 0(a) improved. As a result of these 
technical benefits we have been able to control the disconnected contributions 
and provide statistically significant values for cTttAt = 'mi{N\uu + dd\N) and 
ctq = mi{N\uu + dd- 2ss|A^). 

In the case of the scalar charm content of the nucleon, our statistics was not 
sufficiently large to yield a signal above the statistical noise. We could thus 
only give the bound |(A^|Oc|iV)| < | (iV|0,|A^) |. 

We remark that for phenomenological applications, the relevant quantities are 
actually the RGI quantities ms{N\Os\N) and mc{N\Oc\N) . Given our current 
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statistical accuracy, it is therefore unclear at this moment whether the large 
Yukawa coupling of the charm quark can compensate the smallness of the 
matrix element. 

The most important achievement of this paper is the rather accurate eval- 
uation of the ratio y]^ = 2{N\ss\N) / {N\uu + dd\N) between the strange 
and the light quark content of the nucleon (see Eq. ([3])), for which we find 
Un = 0.082(16)(2). The value we obtain is small compared to estimates from 
chiral perturbation theory but is in line with recent lattice results obtained 
by other groups [2I1I2S1EZII1S] • 

Naturally the results presented in this paper need to be further scrutinized. In 
particular a careful study of the unwanted excited state contamination must 
be carried out to reduce the magnitude of the systematic errors associated to 
these effects. Finally data points at various lattice spacings and pion masses 
are necessary to be able to safely perform an extrapolation to the continuum 
limit and to the physical pion mass. 
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A Scalar density renormalization with maximally twisted Wilson 
quarks 



In this Appendix we want to discuss the renormalization properties of scalar 
quark operators. We will separately discuss the unitary and the mixed action 
case in the setting offered by maximally twisted Wilson fermions [5U] . 

We work here in the so-called physical quark basis and adopt the notations of 
refs. [3TI433] . We recall however that, as it is customary (see e.g. ref. [49]), the 
operator renormalization constants (RC's), being independent of the twisting 
angle (in all mass-independent schemes), are named after the form operators 
take in the twisted quark basis. 



A.l Unitary degenerate doublet 



It has been proved in refs. [30|[3T| that at maximal twistO] in the case of a 
degenerate u, d doublet, with /i„ = — /i^ = /i;, quark mass and scalar density 
renormalize according to the formulae 



,R 



[uu + dd)'^ = 



uu + dd 



I 2 2 2\/^'-n 



(A.l) 
(A.2) 



where Eq. flA.2p is written in the physical quark basis (with r„ = — as a 
consequence of the above chosen values of /x„ and [li) and the last term in 
its r.h.s. represents the mixing of the quark scalar density operator with the 
identity. 

In this paper this term is of no importance because it will be automatically 
subtracted out in the computation of the nucleon matrix elements as described 
in the main text. For this reason, in order not to overload the forthcoming 
formulae, this mixing will not be indicated anymore. 



A.2 Unitary non- degenerate doublet 



For a pair of maximally twisted mass non-degenerate quarks, which for con- 
creteness we name s and c, one finds [311 [33] the more complicated relations 

^ We recall that at maximal twist the flavour group of the lattice theory is the 
so-called 5'C/(2)obiique group with generators {Q\,Q\,Qy}- 
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/if 
/if 



+ e-^ = Zp^rh + Zgh . 



— R 

m 



-1^ 



Zsh 



(A.3) 
(A.4) 



{cc)^ = ^{cc + ss) + ^{cc-ss), (A.5) 
{ss)^ = —^{uu + ss) 2~^^^ ~ ' (^-6) 

where the bare mass parameters rh and e coincide, respectively, with the av- 
erage mass and the mass difference of c and s quarks in free theory. 



A.3 OS valence fermions 



The mass RC of each OS valence fermion in the action is = Zp^. This 
follows from the proof provided in ref. [33] or from the extension of an old 
argument given in ref. [50] that we reproduce in sect. [Bjfor completeness. We 
thus get 

/xg5 = ZpVo5. (A.7) 
We recall that Zp is an even function of the r- Wilson parameter. 

In the philosophy of the "mixed action" approach proposed in ref. [33], a pair 
of mass degenerate OS fermions, denoted (in the so-called physical basis) as 
s_|_ and s_, with opposite values of the Wilson parameter (r^^ = — = 1) 
is introduced to represent the valence s quark, with the understanding that 
no Wick contractions between the fermion s+ and the fermion s_ is allowed. 
This is done also in sect. \2.2\ with /i^ > denoting the bare s quark mass. 

Consider a correlator where besides the strange quark scalar density only 
(renormalized) operators containing no strange quark are present. Then the 
insertion of the renormalized combination (we recall that the divergent mixing 
with the identity must be subtracted out) 

(ss)^ = ^[s+s+ + s„s_], r,^ = -r,_ = l. (A.8) 

is finite, i.e. no new divergences are introduced. The reason for this fact can 
be traced back to the cancellation of chiral violating effects (coming from 
"quark disconnected" - i.e. OZI [SIHS3] violating - diagrams) between the two 
self-contractions ( "loops" ) of the two valence quarks regularized with opposite 
values of r. Alternatively this result can be ascribed to the fact that, having the 
members of the s+, s_ pair opposite values of the r parameters they look like 
a mass degenerate (valence) flavour doublet (in the "physical" quark basis), 
e.g. just as the mass degenerate u and d pair discussed above. Naturally it 
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remains the fact that the theory is not unitary since valence and sea quarks 
are regularized differently. This lack of unitarity leaves behind only O(a^) 
effects [33] • 



B Mass renormalization for OS valence fermions 



For completeness in this appendix we want to explicitly prove the relation f lA.7P 
along the lines of ref. jSU] in a setting where N^>2 valence OS fermions are in- 
troduced over an SU {Nf = 2) (or SU {Nf = 2 + 1 + 1)) maximally twisted sea. 
The renormalization condition (1A.7P is valid for anyone of the A'^^ OS fermions, 
being Zp the renormalization constant of the non-singlet pseudo-scalar quark 
density. 

The key observation to prove Eq. ( 1A.7P is to recall that in order to be entitled 
to use the techniques that are usually employed to derive WTIs, it is necessary 
to have a fully local formulation of the theory. This means that for the purpose 
of dealing with a mixed action case, one has to keep in mind that for each 
OS valence quark in the action a corresponding ghost with equal mass (and 
opposite statistics) has to be introduced. Without loss of generality, for the 
purpose of proving Eq. (IA.7|1 . we can assume that all valence quarks (and 
ghosts) have the same bare (twisted) mass, fios- The generalization to the 
non-degenerate valence quarks is straightforward. 



WTIs for OS fermions 



Let A"^, a = 1,2,...,N^ - 1 be the (non-singlet) axial vector current con- 
structed in terms of only valence fermions (and no ghosts) and let us assume 
that we are exactly at maximal twist. For convenience we shall work in the 
twisted fermion basis where the valence quark mass term has the expression 



LZass = a^Y^f^os Y.iXf{x)il5Xf{x) + ghosts) 



(B.l) 



The argument can be split into four parts 

I - The (bare/lattice) axial WTI between on-shell states reads [51] 

{a\V,A^^iom = 2/ios(«|^'^(0)|/3) - («|X^(0)|/3) , (B.2) 

where 

= xr^X (B.3) 
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and X" is the chiral variation of the OS Wilson term, the exphcit expression 
of which we do not need in this discussion. The only thing we need to know 
about is its mixing pattern (dictated by dimensional argument and the 
symmetry x {nos ~l^os)i where is the formal parity acting on the 
twisted fields, see ref. [30]) which reads 

= {Za - l)V,Al{x) + 2r^x^^osS''{x) + X'^(x) , (B.4) 

with Tjx a (finite) function of g^. Inserting Eq. (lB.4p in (lB.2p . one gets 

(a|Z^V^A;(0)|/3) = 2(1 - Vx)^ios{oi\S%m) " («l^"(0)|/3) , (B.5) 

Since between on-shell states X"- can only give rise to 0(a) terms, the renor- 
malized (continuum looking) WTI 

{a\ZAV,AimP) = 2^gc,(a|5S(0)|/3) + 0(a) . (B.6) 
is immediately obtained by setting 

= ZsS" , (B.7) 
^^os = (1 - Vx)Zs^fios = Z^^fios ■ (B.8) 

Our aim is to prove the relation 

Zp' = il-vx)Zs' (B.9) 

from which the formula 

Zp^ = Z^ (B.IO) 

follows. 

II - To this end we need to extend the previous equations to the case where the 
divergence of the axial current is inserted together with the singlet pseudo- 
scalar density 

^0 = J2(Xfil5Xf + ghosts) (B.ll) 
/ 

which results from considering contributions coming from valence quarks as 
well as the associated ghosts. We note that as the axial current we are con- 
sidering is only made up of valence quarks, it cannot rotate the ghost fields. 
We thus get for the WTI where the operator Pq is inserted 



{a\V,A;{x)P\ym = 2(«|^»(x)|/3)5(x - y) 
+2f,os{a\S%x)P'{ym - (a|X"(x)P°(y)|/3) 



(B.12) 



with external states such that one does not get identically vanishing matrix 
elements. Introducing the decomposition (lB.4p . we first rewrite the previous 
equation in the form 
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{a\ZAV^Al{x)P,{ym = 2(a|^'^(x)|/3)5(x - y) + 

+2(1 - 1lx)^Jios{a\S\x)Po{ym - {a\X\x)P^{ym . (B.13) 

We remark that, when X'^ is inserted with a local operator, it gives rise to 
localized terms plus genuinely 0(a) contributions. Simple symmetry consider- 
ations and the fact that in X" only valence fermions (and not ghosts) appear, 
imply 

{a\X\x)Po{ym = 2cx(a|5"(x)|/3)5(x - y) + 0(a) , (B.14) 

where cx is again a finite function of Qq. Substituting into Eq. (IB.13P and 
neglecting irrelevant (for the present argument) 0(a) terms gives 



{a\ZAV,Al{x)P^\y)\l3) = 2(1 - cx)%(a|^'^(x)|/3)5(x - y) + 

+2^i^s{ol\S^R{^)P^{ym ^ (B.15) 

where we have multiplied both members by Zp^ and used Eqs. flB.7p and fIB.Sp . 
It must be stressed that consistency with continuum WTIs (universality) at 
vanishing iiqs [MIES] requires the identification 

Zpo{l-cx) = Zs. (B.16) 



III - The third step of this analysis is inspired by the discussion carried out 
at the end of sect. 2 of ref. [50]. One notices that by summing over x in both 
members of Eq. flB.15p . one gets at iiqs 7^ the identity 

= 2{a\S'^\ym + 2/.g^(a| ^ S^\x)P^{ym . (B.17) 

X 

as the integral of a divergence vanishes. Once also summed over Eq. f lB.17p 
can be usefully compared to the formula one gets from the obvious identity 

= T^{a\Y.ZAV,Al{xm = T^{a\Y.2f^BsS^\xm , (B.18) 

in which the second equality follows from Eq. (IB. 61) . Indeed by explicitly per- 
forming the derivative with respect to fios, one finds 

= 2(1 - Vx)Zs'{(^\ E S'^'^ixm + 2/ig^(a| ^ S'^'^ix) ^ Po(z/)|/3) • (B.19) 

X X y 

Multiplying this equation by Zp^ and comparing with ( IB. 171) summed over y, 
one gets 

Zp^il-rjx) = Zs. (B.20) 
We remark that Eq. f lB.20p taken together with Eq. ( ]B.16|) entails the some- 
what surprising equality rjx = cx- 
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IV - As the final step of tliis long argument we want to now show that 



% = Zp (B.21) 

which would finally prove Eq. fIB.Qp . The reason for the validity of Eq. flB.2ip is 
that, when one considers the diagrams contributing to Zp^, one realizes that 
the OZI-violating diagrams where a valence quark is self-contracted (closed 
into a "loop" ) is exactly cancelled by the contribution where the corresponding 
ghost, present in Pq? is closed into a "loop". One is thus only left with diagrams 
in which neither valence nor ghost self-contractions appear, hence exactly with 
the diagrams that contribute to the non-singlet pseudo-scalar density RC, Zp, 
where such self-contractions are forbidden by flavour conservation. 
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